Spiec Easi: Sparse InversE Covariance estimation for Ecological Association and Statistical Inference

DUBii 2021

Chargement des librairies

Chargement des données brutes

Transformation

• Transformation log2 centré et mise à l’échelle interquartile pour CPMG et NOESY
  
• Transformation CLR pour OTU (log et centré = center log ratio)
  cf code html métabo

Chargement et filtre des données

Garder uniquement les OTUs présents dans >30% des échantillons (on ne fera rien des autres pour le réseau).
compter nbre de 0 dans chaque colonne et ne retenir que ceux présents dans au moins 30% des échantillons soit 253 sur 846

load("Data_Biom.Rdata")
load("MultiDataSet.Rdata")



otu_trsf <- as.data.frame(otu_table(filter_biom))
count <- apply(otu_trsf, 2, function(x){length(which(x!=0))} )
nb_count <- which(count>=253)
otu_trsf2 <- otu_trsf[,nb_count]
#as.matrix(as.integer(
#otu_table(filter_biom, taxa_are_rows = FALSE) <- otu_trsf2
otu <- as.matrix(otu_trsf2)
taxa <- tax_table(filter_biom)
Sample_ID <- sample_data(filter_biom)

Fil_Biom <- phyloseq(otu_table(otu, taxa_are_rows = FALSE),
  tax_table(taxa),
  sample_data(Sample_ID))
taxa <- as.data.frame(tax_table(Fil_Biom))

filter_biom <- Fil_Biom

Sortie : >huge::huge.roc(se\(est\)path, graph, verbose=FALSE)

True Postive Rate: from 0 to 0.988806 False Positive Rate: from 0 to 0.7952971 Area under Curve: 0.7800105 Maximum F1 Score: 0.9760311 > stars.pr(getOptMerge(se), graph, verbose=FALSE)

True Postive Rate: from 0 to 0.9626866 False Positive Rate: from 0 to 0.1032948 Area under Curve: 0.9778643 Maximum F1 Score: 0.9642207

chargement des simulations effectuées avec différentes méthodes

MB, Glasso et SparCC
Exemple de commande :
se.mb.biom_1 <- spiec.easi(otu_table(filter_biom), method=‘mb’, lambda.min.ratio=1e-3, nlambda=20, pulsar.params=list(rep.num=50, ncores=16))
se.gl.biom_1 <- spiec.easi(otu_table(filter_biom), method=‘glasso’, lambda.min.ratio=1e-3, nlambda=20,pulsar.params=list(rep.num=50, ncores=16))
réalisé sur différents lambda.min.ratio
λ=1e-3 λ=1e-2 λ=0.05 λ=0.1 λ=0.15 λ=0.2 λ=0.3

#load("spiec_easi.Rdata")
#load("spiec_easi2.Rdata")
#load("spiec_easi_filt.Rdata")
#se.mb.biom_init <- se.mb.biom
#se.gl.biom_init <-  se.gl.biom
#load("spiec_easi_filt_teste1.Rdata")
#load("spiec_easi_filt_teste2.Rdata")
#load("spiec_easi_p1.Rdata")


load("spiec_easi_gl_mb.Rdata")

## Define arbitrary threshold for SparCC correlation matrix for the graph
sparcc.graph <- abs(sparcc.biom$Cor) >= 0.3
diag(sparcc.graph) <- 0

sparcc.graph <- Matrix(sparcc.graph, sparse=TRUE)
## Create igraph objects

#lambda.min.ratio=1e-3  
ig.mb_1     <- adj2igraph(getRefit(se.mb.biom_1))
ig.gl_1     <- adj2igraph(getRefit(se.gl.biom_1))
ig.sparcc_1 <- adj2igraph(sparcc.graph)

#lambda.min.ratio=1e-2 val par défaut
ig.mb     <- adj2igraph(getRefit(se.mb.biom))
ig.gl     <- adj2igraph(getRefit(se.gl.biom))
ig.sparcc <- adj2igraph(sparcc.graph)

#lambda.min.ratio=5e-2
ig.mb005     <- adj2igraph(getRefit(se.mb.biom005))
ig.gl005     <- adj2igraph(getRefit(se.gl.biom005))
ig.sparcc005 <- adj2igraph(sparcc.graph)

#lambda.min.ratio=0.1
ig.mb01     <- adj2igraph(getRefit(se.mb.biom01))
ig.gl01     <- adj2igraph(getRefit(se.gl.biom01))
ig.sparcc01 <- adj2igraph(sparcc.graph)

#lambda.min.ratio=0.15
ig.mb015    <- adj2igraph(getRefit(se.mb.biom015))
ig.gl015    <- adj2igraph(getRefit(se.gl.biom015))
ig.sparcc015 <- adj2igraph(sparcc.graph)

#lambda.min.ratio=0.2
ig.mb02     <- adj2igraph(getRefit(se.mb.biom02))
ig.gl02     <- adj2igraph(getRefit(se.gl.biom02))
ig.sparcc02 <- adj2igraph(sparcc.graph)


#lambda.min.ratio=0.3
ig.mb03     <- adj2igraph(getRefit(se.mb.biom03))
ig.gl03    <- adj2igraph(getRefit(se.gl.biom03))
ig.sparcc03 <- adj2igraph(sparcc.graph)

Visualisation

## set size of vertex proportional to clr-mean
vsize    <- rowMeans(clr(otu_table(filter_biom), 1))+6
biom.coord_1 <- layout.fruchterman.reingold(ig.mb_1)
biom.coord <- layout.fruchterman.reingold(ig.mb)
biom.coord005 <- layout.fruchterman.reingold(ig.mb005)
biom.coord01 <- layout.fruchterman.reingold(ig.mb01)
biom.coord015 <- layout.fruchterman.reingold(ig.mb015)
biom.coord02 <- layout.fruchterman.reingold(ig.mb02)
biom.coord03 <- layout.fruchterman.reingold(ig.mb03)

biom.coord_g1 <- layout.fruchterman.reingold(ig.gl_1)
biom.coordg <- layout.fruchterman.reingold(ig.gl)
biom.coordg005 <- layout.fruchterman.reingold(ig.gl005)
biom.coordg01 <- layout.fruchterman.reingold(ig.gl01)
biom.coordg015 <- layout.fruchterman.reingold(ig.gl015)
biom.coordg02 <- layout.fruchterman.reingold(ig.gl02)
biom.coordg03 <- layout.fruchterman.reingold(ig.gl03)


sparcc.graph1 <- abs(sparcc.biom$Cor) >= 0.01
diag(sparcc.graph) <- 0
sparcc.graph1 <- Matrix(sparcc.graph, sparse=TRUE)
ig.sparcct1 <- adj2igraph(sparcc.graph1)

sparcc.graph2 <- abs(sparcc.biom$Cor) >= 0.9
diag(sparcc.graph) <- 0
sparcc.graph2 <- Matrix(sparcc.graph, sparse=TRUE)
ig.sparcct2 <- adj2igraph(sparcc.graph2)

par(mfrow=c(1,3))
plot(ig.sparcc, layout=biom.coord, vertex.size=vsize, vertex.label=NA, main="sparcc 0.3 layoutmb")
plot(ig.sparcct1, layout=biom.coord, vertex.size=vsize, vertex.label=NA, main="sparcc 0.01")
plot(ig.sparcct2, layout=biom.coord, vertex.size=vsize, vertex.label=NA, main="sparcc 0.9")

par(mfrow=c(1,3))
plot(ig.sparcc, layout=biom.coordg, vertex.size=vsize, vertex.label=NA, main="sparcc 0.3 layoutglasso")
plot(ig.sparcct1, layout=biom.coordg, vertex.size=vsize, vertex.label=NA, main="sparcc 0.01")
plot(ig.sparcct2, layout=biom.coordg, vertex.size=vsize, vertex.label=NA, main="sparcc 0.9")

par(mfrow=c(1,3))
plot(ig.mb_1, layout=biom.coord_1, vertex.size=vsize, vertex.label=NA, main="MB λ=1e-3  ")
plot(ig.gl_1, layout=biom.coord_g1, vertex.size=vsize, vertex.label=NA, main="glasso")
plot(ig.sparcc_1, layout=biom.coord_1, vertex.size=vsize, vertex.label=NA, main="sparcc")

par(mfrow=c(1,3))
plot(ig.mb, layout=biom.coord, vertex.size=vsize, vertex.label=NA, main="MB λ=1e-2  ")
plot(ig.gl, layout=biom.coordg, vertex.size=vsize, vertex.label=NA, main="glasso")
plot(ig.sparcc, layout=biom.coord, vertex.size=vsize, vertex.label=NA, main="sparcc")

par(mfrow=c(1,3))
plot(ig.mb005, layout=biom.coord005, vertex.size=vsize, vertex.label=NA, main="MB λ=0.05    ")
plot(ig.gl005, layout=biom.coordg005, vertex.size=vsize, vertex.label=NA, main="glasso")
plot(ig.sparcc005, layout=biom.coord005, vertex.size=vsize, vertex.label=NA, main="sparcc")

par(mfrow=c(1,3))
plot(ig.mb01, layout=biom.coord01, vertex.size=vsize, vertex.label=NA, main="MB λ=0.1   ")
plot(ig.gl01, layout=biom.coordg01, vertex.size=vsize, vertex.label=NA, main="glasso")
plot(ig.sparcc01, layout=biom.coord01, vertex.size=vsize, vertex.label=NA, main="sparcc")

par(mfrow=c(1,3))
plot(ig.mb015, layout=biom.coord015, vertex.size=vsize, vertex.label=NA, main="MB λ=0.15    ")
plot(ig.gl015, layout=biom.coordg015, vertex.size=vsize, vertex.label=NA, main="glasso")
plot(ig.sparcc, layout=biom.coord015, vertex.size=vsize, vertex.label=NA, main="sparcc")

par(mfrow=c(1,3))
plot(ig.mb02, layout=biom.coord02, vertex.size=vsize, vertex.label=NA, main="MB λ=0.2   ")
plot(ig.gl02, layout=biom.coordg02, vertex.size=vsize, vertex.label=NA, main="glasso")
plot(ig.sparcc02, layout=biom.coord02, vertex.size=vsize, vertex.label=NA, main="sparcc")

par(mfrow=c(1,3))
plot(ig.mb03, layout=biom.coord03, vertex.size=vsize, vertex.label=NA, main="MB λ=0.3   ")
plot(ig.gl03, layout=biom.coordg03, vertex.size=vsize, vertex.label=NA, main="glasso")
plot(ig.sparcc03, layout=biom.coord03, vertex.size=vsize, vertex.label=NA, main="sparcc")

Evaluation du poids des liaisons du reseau (edge)

secor_1  <- cov2cor(getOptCov(se.gl.biom_1))
sebeta_1 <- symBeta(getOptBeta(se.mb.biom_1), mode='maxabs')
elist.gl_1     <- summary(triu(secor_1*getRefit(se.gl.biom_1), k=1))
elist.mb_1     <- summary(sebeta_1)
elist.sparcc <- summary(sparcc.graph*sparcc.biom$Cor)
hist(elist.mb_1[,3], main='λ=1e-3  grey=MB\n green=sparcc, red=Glasso', xlab='edge weights')
hist(elist.sparcc[,3], add=TRUE, col='forestgreen')
hist(elist.gl_1[,3], add=TRUE, col='red')

secor  <- cov2cor(getOptCov(se.gl.biom))
sebeta <- symBeta(getOptBeta(se.mb.biom), mode='maxabs')
elist.gl     <- summary(triu(secor*getRefit(se.gl.biom), k=1))
elist.mb     <- summary(sebeta)
elist.sparcc <- summary(sparcc.graph*sparcc.biom$Cor)
hist(elist.mb[,3], main='λ=1e-2  grey=MB\n green=sparcc, red=Glasso', xlab='edge weights')
hist(elist.sparcc[,3], add=TRUE, col='forestgreen')
hist(elist.gl[,3], add=TRUE, col='red')

secor005  <- cov2cor(getOptCov(se.gl.biom005))
sebeta005 <- symBeta(getOptBeta(se.mb.biom005), mode='maxabs')
elist.gl005     <- summary(triu(secor005*getRefit(se.gl.biom005), k=1))
elist.mb005     <- summary(sebeta005)
elist.sparcc <- summary(sparcc.graph*sparcc.biom$Cor)
hist(elist.mb005[,3],ylim=c(1,400), main='λ= 0.05   grey=MB\n green=sparcc, red=Glasso', xlab='edge weights')
hist(elist.sparcc[,3], add=TRUE, col='forestgreen')
hist(elist.gl005[,3], add=TRUE, col='red')

secor01  <- cov2cor(getOptCov(se.gl.biom01))
sebeta01 <- symBeta(getOptBeta(se.mb.biom01), mode='maxabs')
elist.gl01     <- summary(triu(secor005*getRefit(se.gl.biom01), k=1))
elist.mb01     <- summary(sebeta01)
elist.sparcc <- summary(sparcc.graph*sparcc.biom$Cor)
hist(elist.mb01[,3], main='λ=0.1 grey=MB\n green=sparcc, red=Glasso', xlab='edge weights')
hist(elist.sparcc[,3], add=TRUE, col='forestgreen')
hist(elist.gl01[,3], add=TRUE, col='red')

secor015  <- cov2cor(getOptCov(se.gl.biom015))
sebeta015 <- symBeta(getOptBeta(se.mb.biom015), mode='maxabs')
elist.gl015     <- summary(triu(secor015*getRefit(se.gl.biom015), k=1))
elist.mb015     <- summary(sebeta015)
elist.sparcc <- summary(sparcc.graph*sparcc.biom$Cor)
hist(elist.mb015[,3], main='λ=0.15 grey=MB\n green=sparcc, red=Glasso', xlab='edge weights')
hist(elist.sparcc[,3], add=TRUE, col='forestgreen')
hist(elist.gl015[,3], add=TRUE, col='red')

secor02  <- cov2cor(getOptCov(se.gl.biom02))
sebeta02 <- symBeta(getOptBeta(se.mb.biom02), mode='maxabs')
elist.gl02     <- summary(triu(secor02*getRefit(se.gl.biom02), k=1))
elist.mb02     <- summary(sebeta02)
elist.sparcc <- summary(sparcc.graph*sparcc.biom$Cor)
hist(elist.mb02[,3], main='λ=0.2 grey=MB\n green=sparcc, red=Glasso', xlab='edge weights')
hist(elist.sparcc[,3], add=TRUE, col='forestgreen')
hist(elist.gl02[,3], add=TRUE, col='red')

secor03  <- cov2cor(getOptCov(se.gl.biom03))
sebeta03 <- symBeta(getOptBeta(se.mb.biom03), mode='maxabs')
elist.gl03     <- summary(triu(secor*getRefit(se.gl.biom03), k=1))
elist.mb03     <- summary(sebeta03)
elist.sparcc <- summary(sparcc.graph*sparcc.biom$Cor)
hist(elist.mb03[,3], main='λ=0.3  grey=MB\n green=sparcc, red=Glasso', xlab='edge weights')
hist(elist.sparcc[,3], add=TRUE, col='forestgreen')
hist(elist.gl03[,3], add=TRUE, col='red')

degrés de statistique des réseaux inférrés par quelques méthodes

dd.gl01     <- degree.distribution(ig.gl01)
dd.mb01     <- degree.distribution(ig.mb01)
dd.sparcc <- degree.distribution(ig.sparcc)
plot(0:(length(dd.sparcc)-1), dd.sparcc, ylim=c(0,.35), type='b',
      ylab="Frequency", xlab="Degree", main="Degree Distributions λ=0.1 ")
points(0:(length(dd.gl01)-1), dd.gl01, col="red" , type='b')
points(0:(length(dd.mb01)-1), dd.mb01, col="forestgreen", type='b')
legend("topright", c("MB", "glasso", "sparcc"),
        col=c("forestgreen", "red", "black"), pch=1, lty=1)

dd.gl015     <- degree.distribution(ig.gl015)
dd.mb015     <- degree.distribution(ig.mb015)
dd.sparcc <- degree.distribution(ig.sparcc)
plot(0:(length(dd.sparcc)-1), dd.sparcc, ylim=c(0,.35), type='b',
      ylab="Frequency", xlab="Degree", main="Degree Distributions λ=0.15")
points(0:(length(dd.gl015)-1), dd.gl015, col="red" , type='b')
points(0:(length(dd.mb015)-1), dd.mb015, col="forestgreen", type='b')
legend("topright", c("MB", "glasso", "sparcc"),
        col=c("forestgreen", "red", "black"), pch=1, lty=1)

dd.gl02     <- degree.distribution(ig.gl02)
dd.mb02     <- degree.distribution(ig.mb02)
dd.sparcc <- degree.distribution(ig.sparcc)
plot(0:(length(dd.sparcc)-1), dd.sparcc, ylim=c(0,.35), type='b',
      ylab="Frequency", xlab="Degree", main="Degree Distributions λ=0.2")
points(0:(length(dd.gl02)-1), dd.gl02, col="red" , type='b')
points(0:(length(dd.mb02)-1), dd.mb02, col="forestgreen", type='b')
legend("topright", c("mb", "glasso", "sparcc"),
        col=c("forestgreen", "red", "black"), pch=1, lty=1)

Visualisation du reseau par famille

Methode MB

## Load round 2 of American gut project
# data('amgut2.filt.phy')
# table_ex <- as.data.frame(tax_table(amgut2.filt.phy))
# se.mb.amgut2 <- spiec.easi(filter_biom, method='mb', lambda.min.ratio=1e-2,
#                            nlambda=20, pulsar.params=list(rep.num=50))
# table_otu <- as.data.frame(tax_table(filter_biom))

#load("param_spiec_clut1.Rdata")
#load("param_spiec_filt2.Rdata")
#lambda.min.ratio=0.1

biom.coord_1 <- layout.fruchterman.reingold(ig.mb_1)
biom.coord <- layout.fruchterman.reingold(ig.mb)
biom.coord005 <- layout.fruchterman.reingold(ig.mb005)
biom.coord01 <- layout.fruchterman.reingold(ig.mb01)
biom.coord015 <- layout.fruchterman.reingold(ig.mb015)
biom.coord02 <- layout.fruchterman.reingold(ig.mb02)
biom.coord03 <- layout.fruchterman.reingold(ig.mb03)

ig_1.mb_1 <- adj2igraph(getRefit(se.mb.biom_1), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.mb_1, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=1e-3" )

ig_1.mb <- adj2igraph(getRefit(se.mb.biom), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.mb, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=1e-2")

ig_1.mb005 <- adj2igraph(getRefit(se.mb.biom005), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.mb005, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=0.005")

ig_1.mb01 <- adj2igraph(getRefit(se.mb.biom01), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.mb01, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=0.1")

ig_1.mb015 <- adj2igraph(getRefit(se.mb.biom015), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.mb015, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=0.15")

ig_1.mb02 <- adj2igraph(getRefit(se.mb.biom02), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.mb01, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=0.2")

ig_1.mb03 <- adj2igraph(getRefit(se.mb.biom03), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.mb03, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=0.3")

# library(phyloseq)
# ## Load round 2 of American gut project
# data('amgut2.filt.phy')
# se.mb.amgut2 <- spiec.easi(amgut2.filt.phy, method='mb', lambda.min.ratio=1e-2,
#                            nlambda=20, pulsar.params=list(rep.num=50))
# ig2.mb <- adj2igraph(getRefit(se.mb.amgut2),  vertex.attr=list(name=taxa_names(amgut2.filt.phy)))
# plot_network(ig2.mb, amgut2.filt.phy, type='taxa', color="Rank3")

Méthode Glasso

biom.coord_1 <- layout.fruchterman.reingold(ig.gl_1)
biom.coord <- layout.fruchterman.reingold(ig.gl)
biom.coord005 <- layout.fruchterman.reingold(ig.gl005)
biom.coord01 <- layout.fruchterman.reingold(ig.gl01)
biom.coord015 <- layout.fruchterman.reingold(ig.gl015)
biom.coord02 <- layout.fruchterman.reingold(ig.gl02)
biom.coord03 <- layout.fruchterman.reingold(ig.gl03)


ig_1.gl_1 <- adj2igraph(getRefit(se.gl.biom_1), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.gl_1, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=1e-3" )

ig_1.gl <- adj2igraph(getRefit(se.gl.biom), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.gl, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=1e-2")

ig_1.gl005 <- adj2igraph(getRefit(se.gl.biom005), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.gl005, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=0.05")

ig_1.gl01 <- adj2igraph(getRefit(se.gl.biom01), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.gl01, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=0.1")

ig_1.gl015 <- adj2igraph(getRefit(se.gl.biom015), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.gl015, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=0.15")

ig_1.gl02 <- adj2igraph(getRefit(se.gl.biom02), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.gl01, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=0.2")

ig_1.gl03 <- adj2igraph(getRefit(se.gl.biom03), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig_1.gl03, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="λ=0.3")

test à l’espèce

Trop d’espèce impossible de visualiser le reseau. Il faudrait ne garder que les 10 plus representer

plot_network(ig_1.mb, filter_biom, type='taxa', color="Species", label = NULL, line_weight=0.4, title ="λ=1e-2")

Learning latent variable graphical models

load("param_spiec_slr.Rdata")

#se.slr10.biom <- spiec.easi(otu_table(filter_biom), method='slr', r=10, lambda.min.ratio=1e-2,
#                           nlambda=20, pulsar.params=list(rep.num=50, ncores=16))

#Différents param testés slr 10 et 20 lambda e-2 et 0.1



ig.slr10    <- adj2igraph(getRefit(se.slr10.biom))
biom.coordslr10 <- layout.fruchterman.reingold(ig.slr10)

ig.slr20    <- adj2igraph(getRefit(se1.slr20.biom))
biom.coordslr20 <- layout.fruchterman.reingold(ig.slr20)

ig.slr101    <- adj2igraph(getRefit(se.slr10.biom1))
biom.coordslr101 <- layout.fruchterman.reingold(ig.slr101)

ig.slr201    <- adj2igraph(getRefit(se1.slr20.biom1))
biom.coordslr201 <- layout.fruchterman.reingold(ig.slr201)

vsize    <- rowMeans(clr(otu_table(filter_biom), 1))+6


par(mfrow=c(2,2))
plot(ig.slr10, layout=biom.coordslr10, vertex.size=vsize, vertex.label=NA, main="slr 10 λ=1e-2")
plot(ig.slr20, layout=biom.coordslr20, vertex.size=vsize, vertex.label=NA, main="slr 20 λ=1e-2")
plot(ig.slr101, layout=biom.coordslr101, vertex.size=vsize, vertex.label=NA, main="slr 10 λ=0.1")
plot(ig.slr101, layout=biom.coordslr201, vertex.size=vsize, vertex.label=NA, main="slr 20 λ=0.1")

ig1.slr10 <- adj2igraph(getRefit(se.slr10.biom), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig1.slr10, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="slr 10 λ=1e-2" )

ig1.slr20 <- adj2igraph(getRefit(se1.slr20.biom), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig1.slr20, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="slr 20 λ=1e-2" )

ig1.slr101 <- adj2igraph(getRefit(se.slr10.biom1), vertex.attr=list(name=taxa_names(filter_biom)))
plot_network(ig1.slr101, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="slr 10 λ=0.1" )

# ig1.slr201 <- adj2igraph(getRefit(se1.slr20.biom1), vertex.attr=list(name=taxa_names(filter_biom)))
# plot_network(ig1.slr201, filter_biom, type='taxa', color="Phylum", label = NULL, line_weight=0.4, title ="slr 20 λ=0.1" )